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Abstract 

It  is  known  that  the  exact  analytic  solutions  of  wave  scattering  by  a  circular  cylin¬ 
der,  when  they  exist,  are  not  in  a  closed  form  but  in  infinite  series  which  converges 
slowly  for  high  frequency  waves.  In  this  paper,  we  present  a  fast  numerical  solution 
for  the  scattering  problem  in  which  the  Boundary  Integral  Equations,  reformulated 
from  the  Helmholtz  equation,  are  solved  using  a  Fourier  spectral  method.  It  is  shown 
that  the  special  geometry  considered  here  allows  the  implementation  of  the  spectral 
method  to  be  simple  and  very  efficient.  The  present  method  differs  from  previous  ap¬ 
proaches  in  that  the  singularities  of  the  integral  kernels  are  removed  and  dealt  with 
accurately.  The  proposed  method  preserves  the  spectral  accuracy  and  is  shown  to  have 
an  exponential  rate  of  convergence.  Aspects  of  efficier*  implementation  using  FFT  are 
discussed.  Moreover,  the  boundary  integral  equations  of  combined  single  and  double¬ 
layer  representation  are  used  in  the  present  paper.  This  ensures  the  uniqueness  of  the 
numerical  solution  for  the  scattering  problem  at  all  frequencies.  Although  a  strongly 
singular  kernel  is  encountered  for  the  Neumann  boundary  conditions,  we  show  that  the 
hypersingularity  can  be  handled  easily  in  the  spectral  method.  Numerical  examples 
that  demonstrate  the  validity  of  the  method  are  also  presented. 
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NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Science  and 
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INTRODUCTION 


The  exact  analytic  solutions  of  wave  scattering  by  a  circular  cylinder,  obtainable 
for  simple  incident  waves,  are  not  in  a  closed  form  but  in  infinite  series  of  Bessel  and 
Hankel  functions  of  increasing  orders.  Such  solutions  converge  slowly,  especially  for 
high  frequency  waves,  which  render  their  numerical  evaluation  ineflSicient.  This  paper 
presents  a  fast  numerical  solution  of  wave  scattering  that  only  requires  the  computation 
of  Bessel  and  Hankel  functions  of  order  zero.  Furthermore,  the  numerical  solution  is 
valid  for  any  form  of  the  incident  waves  of  all  frequencies. 

When  developing  numerical  solutions,  wave  scattering  problems  are  often  conve¬ 
niently  formulated  in  Boundary  Integral  Equations  (BIE)*.  The  advantages  of  the 
Boundary  Integral  Equation  Method  (BIEM)  include  reducing  the  dimension  of  the 
problem  and  transforming  an  infinite  domain  to  finite  boundaries  in  which  the  far  field 
radiation  condition  is  satisfied  automatically.  The  Boundary  Integral  Equations  are 
commonly  solved  computationally  by  the  Boundary  Element  Methods  (BEM)^.  In 
this  method  the  boundary  is  divided  into  finite  elements  and  integrations  over  each 
boundary  element  are  approximated  by  quadratures,  e.g.  the  linear  elements. 

In  this  paper,  we  develop  a  Spectral  Method  of  solving  the  Boundary  Integral 
Equations,  reformulated  from  the  Helmholtz  equation,  for  numerical  solutions  of  wave 
scattering  by  a  circular  cylinder.  Previously,  for  this  special  geometry,  a  “fast  numerical 
method”  based  on  the  Fourier  approximations  has  been  formulated  by  Bojarski^,  who 


approach  was  more  efficient  than  directly  evaluating  the  infinite  series  of  the  exact 
solutions.  Indeed,  the  exaict  solutions  contain  Bessel  and  Hankel  functions  of  higher 
orders  whose  numerical  evaluation  is  more  difficult  and  costly  as  the  order  increases. 
Recently  a  similar  approach  has  been  used  and  extended  by  Schuster^  for  a  wave 
transmission  problem  of  concentric  cylinders. 

In  the  present  paper,  we  point  out  that  the  numerical  formulations  given  previously 
are  not  achieving  the  optimal  accuracy  of  the  Fourier  spectral  methods.  It  is  known 
that,  although  any  periodic  function  can  be  approximated  by  a  truncated  Fourier  se¬ 
ries,  the  rate  of  convergence  of  such  an  approximation  depends  on  its  smoothness. 
Unfortunately,  the  integral  kernels  for  the  Helmholtz  equation  are  not  smooth.  In  par¬ 
ticular,  the  2-D  Green’s  function  of  the  Helmholtz  equation,  appearing  in  the  integral 
equations,  possesses  a  logarithmic  singularity.  Furthermore,  the  normal  derivative  of 
the  Green’s  function  also  contains  a  term  involving  the  logarithmic  function.  The  non¬ 
smoothness  of  the  integral  kernels,  however,  was  not  explicitly  treated  in  the  previous 
formulations.  It  will  be  seen  that  it  is  critical  to  remove  the  non-smoothness  of  the 
integral  kernels  in  order  to  achieve  fast  convergence  in  the  Fourier  spectral  formulation. 
By  a  proper  treatment  of  the  singularities,  the  present  numerical  formulation  yields 
accurate  solutions  with  significantly  fewer  datum  points.  Moreover,  the  boundary  in¬ 
tegral  equations  of  combined  single  and  double-layer  representation  are  used  in  the 
present  paper.  This  ensures  the  uniqueness  of  the  numerical  solution  for  the  scatter¬ 
ing  problem  at  all  frequencies^’®.  Although  a  combined  layer  formulation  results  in 
a  strongly  singular  kernel  for  the  Neumann  boundary  conditions,  we  show  that  the 
hypersingularity  is  handled  easily  in  the  spectral  method. 

In  the  next  section,  the  formulations  of  the  Boundary  Integral  Equations  for  wave 
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scattering  problems  are  given.  Then,  in  sections  II  and  III,  the  Fourier  spectral  methods 
for  the  Dirichlet  and  Neumann  boundary  conditions  are  presented.  Numerical  results 
are  shown  in  section  IV.  Section  V  contains  the  conclusions.  Some  analytic  results  are 
also  given  in  the  Appendix. 

I.  BOUNDARY  INTEGRAL  EQUATIONS 

Let  us  consider  wave  scattering  by  a  circular  cylinder  F  of  radius  a.  The  wave 
equation  for  the  scattered  function  <f>  with  assumed  time  dependency  of  is  reduced 
to  the  Helmholtz  equation 

V2<^  +  »c2^  =  0  (1) 

where  k  =  ujc  (c  is  the  wave  speed)  and  is  the  2-D  Laplace  operator  = 
.  The  boundary  condition  considered  in  this  paper  will  be  one  of  the  following 

types  : 

Dirichlet{soft)  :  <i>{r)  =  6(r )  On  F 

or 

d4> 

Neumann{hard)  :  ^ 

The  Helmholtz  equation  (1)  together  with  the  boundary  condition  can  be  refor¬ 
mulated  into  a  Boundary  Integral  Equation.  This  can  be  done  in  various  ways^’®.  For 
scattering  problems  considered  in  the  present  paper,  we  use  a  combination  of  single 
and  double-layer  formulation  in  which  the  solution  <i>  at  any  point  r'  in  the  scattered 
field  is  represented  by  an  integral  on  the  boundary  as® 

')  =  /  (^  -  ‘ig)  /(r  )rfr  (2) 


dy'^ 
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where  is  any  real  number  such  that 


rjRe{K)  >  0 


The  use  of  a  combined  formulation  ensures  the  uniqueness  of  the  numerical  solution 
for  exterior  problems*’®.  In  (2),  /(f)  is  an  unknown  layer  distribution  function  and 
the  Green’s  function  G(f,  f '),  whose  form  will  be  given  later,  satisfies  the  following 


equation 


V^G  +  K'^G  =  -6{r-f>) 


Here  the  normal  derivative  is  assumed  to  be  taken  in  the  direction  outward  from 

an 

the  cylinder. 

The  Boundary  Integral  Equation  associated  with  the  layer  representation  (2)  is® 


for  Dirichlet  boundary  conditions  and 

-  •-'9  nrv^ = Mr-;)  («) 

for  Neumann  boundary  conditions,  respectively.  In  (4a)  and  (4b),  fp  denotes  the  bound¬ 
ary  points.  After  the  layer  distribution  function  /  has  been  solved  from  the  integral 
equation  (4a)  or  (4b),  the  solution  of  the  Helmholtz  equation  <j>  is  found  by  the  bound¬ 
ary  integral  (2). 

Now  for  a  circular  cylinder  of  radius  a,  the  boundary  contour  can  be  expressed  as 


fp(d)  =  (acos0,asin^),  0  <  ^  <  27r 


The  normal  vector  to  be  used  in  (4a)  and  (4b)  is  n  =  (cos  0,  sin  0). 
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The  Green’s  function  and  its  normal  derivative  are®’® 

a  =  («|r-r(«)  -  fr(l»')|)  =  j//*"  (2«.  Lin  ^  )  (6) 

and 


dG  _ 

(r-r(«)  -  fr(«'))  ■  n 

dn 

|rp(<?)  -  rp(0')l 

= 

4  (^2«a  |sin  ^  j 

.  0-9' 
sm  — - — 

2 

(7) 

in  which  we  have  used  the  fact  that  Irp(fl)  —  rp{0')|  =  2a  sin  . 

It  is  important  to  note  here  that  G  and  ^  are  functions  of  0  —  6' .  As  will  be  seen 
later,  this  allows  the  implementation  of  the  Fourier  spectral  method  to  take  a  simple 
form. 

We  thus  express  the  boundary  integral  equation  (4a)  for  the  Dirichlet  boundary 
conditions  as 

and  equation  (4b)  for  the  Neumann  boundary  conditions  as 

For  clarity,  the  dependencies  on  9  and  6'  have  been  expressed  explicitly  in  (8a)  and 
(8b). 

In  the  next  two  sections,  we  give  the  numerical  formulations  of  solving  the  inte¬ 
gral  equations  (8a)  and  (8b)  by  a  Fourier  spectral  method.  Since  different  types  of 
singularities  are  encountered,  the  two  equations  will  be  dealt  with  separately. 


il.  SPECTRAL  METHOD  FOR  DIRICHLET  BOUNDARY  CONDITIONS 
A.  Formulation 


Let  the  layer  distribution  function  f{6)  and  the  boundary  condition  b{6)  be  ap¬ 
proximated  by  the  trimcated  Fourier  series  as 

N/2-1 

/(»)=  E  (») 

n=-N/2 

N/2~l 

6(0)=  52  ('») 

n=-N/2 

where  6*  are  obtained  by  the  FFT  from  prescribed  boundary  condition  and  fn  are  the 
unknown  coefficients.  In  (9)  and  (10),  the  particular  form  of  truncated  Fourier  series 
has  been  taken  for  the  convenience  of  applying  FFT  programs. 

Substituting  (9)  and  (10)  into  the  boundary  integral  equation  for  the  Dirichlet 
boundary  conditions  (8a),  we  get 


N/2-l 

5  E  . 

n=-Nl2 


N/2-\ 

/.e‘-'+  E 

n=-N/2 


H"(; 


1(0  -0')-  irjG(e 


-d')^ 


Nl2-\ 

=  S 

n=:-N/2 

(11) 


For  simplicity,  let 


X  =  9  —  ^ 


By  equating  the  coefficients  of  c”*^,  equation  (11)  is  e2isily  reduced  to 

!/«  +  /»  (^(x)  -  WG(x))  e™-o  ix  =  K  (12) 

for  -iV/2  <  n  <  iV/2  -  1. 

It  is  seen  that  the  integral  appearing  in  (12)  are  related  to  the  Fourier  coefficients 
of  ^(x)  amd  G{x).  From  (6)  and  (7),  it  is  also  clear  that  both  are  periodic  functions  of 
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X,  with  a  period  of  2ir.  Thus  if  we  let  G(x)  and  ^(x)  be  approximated  by  truncated 


Fourier  series  as 

N/2-l 

G(x)= 

n=-N/2 

n=-JV/2 

then,  the  integral  in  (12)  equals  to  2ira{hn  —  irign)-  It  follows  that 


"f"  (^n  ~  *V9n)  — 


Therefore,  the  Fourier  coefficients  of  the  layer  distribution  function  f{9)  are  ob¬ 


tained  explicitly  as 


In  -  1 


^  -|-27ra  {hn- irign) 

The  above  equation  shows  that  once  the  Fourier  coefficients  of  G(x)  and  ^{x) 
have  been  found,  the  layer  distribution  function  f{B)  is  known  immediately. 

Actually,  the  Fourier  coefficients  of  G{x)  and  ^(a:)  can  be  found  in  exact  form 
using  higher  order  Bessel  and  Hankel  functions.  They  are  derived  in  Appendix  A. 
Nonetheless,  the  numerical  evaluation  of  the  exact  expressions  becomes  more  ineffective 
and  costly  ais  the  order  of  the  special  functions  increases.  In  what  follows  we  give  the 
numerical  method  that  computes  the  Fourier  coefficients  gn  and  hn  accurately  and 
efficiently. 


B.  Computation  of  gn  and  hn 

In  general,  the  Fourier  coefficients  of  a  periodic  function  can  be  obtained  efficiently 
by  using  a  Fast  Fourier  Transform  algorithm  (FFT).  However,  the  accuracy  of  the 
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Fourier  coefficients  computed  by  the  FFT  using  a  given  number  of  datum  points  de¬ 
pends  on  the  smoothness  of  the  function.  Only  when  the  function  is  infinitely  smooth 
(i.e.  infinitely  differentiable),  the  error  of  Fourier  coefficients  computed  by  FFT  de¬ 
cays  faster  th^m  any  power  of  l/TV,  where  N  is  the  number  of  datum  points.  Such  a 
convergence  is  often  referred  to  as  an  exponential  convergence  and  the  method  is  said 
to  have  spectral  accuracy^’*.  Our  aim  here  is  to  compute  and  h„  by  the  FFT  with 
spectral  accuracy  even  though  the  functions  G  and  ^  are  not  smooth. 

In  the  numerical  approaches  proposed  previously^’^,  the  Fourier  coefficients  gn  and 
hn  were  computed  directly  as  the  FFT  of  the  G{x)  and  ^(x)  respectively.  However, 
the  Green’s  function  G(x)  has  a  logarithmic  singularity  at  a;  =  0,  where  6  =  $',  due 
to  the  Hankel  function  of  order  zero  in  (6),  and  its  Fourier  series  converges  at  the  rate 
of  \/N.  Thus  direct  computation  of  from  G{x)  using  FFT  yields  results  whose 
accuracy  is  only  comparable  to  a  first  order  method.  Furthermore,  the  function  f^(x) 
also  has  a  non-smooth  derivative  at  x  =  0,  and  its  Fourier  series  converges  at  the 
rate  of  1/^^.  Thus  direct  computation  of  h„  from  f^{x)  is  only  comparable  to  a 
third  order  method.  Alternatively,  as  will  be  shown  below,  by  properly  treating  the 
non-smoothness  of  (?(x)  and  ^(x),  g„  and  hn  are  computed  with  spectral  accuracy. 

To  examine  the  singularity  of  G(x),  we  note  that 

G(x)  =  ^^0*^  ^2/ca  sin  ^  [*^0  ^2/ca  jsin  ^2/ca  sin  ^  ^ 

in  which  Jo  and  Vo  are  the  zeroth  order  Bessel  functions  of  the  first  and  second  kind. 


respectively.  Using  the  asymptotic 

,2 


series  lor  small  arguments. 


^»W  =  '-T  +  64-- 

y.W  =  fln(f)j.M+^AW+|i 
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It  follows  that,  for  |x|  small, 


G(x)  =  -^  (««  II)  *^0  (^Ka  sin  ||)  “  ^  +  |  + 

in  which  0{x^)  represents  a  power  series  in  x^,  and  7  is  the  Euler’s  constant,  7  = 
0.577215...  .  To  compute  the  Fourier  coefficients  of  G{x)  efficiently  and  accurately,  we 
note  that  the  Fourier  series  of  the  logarithmic  periodic  function  In  (/ca  |sin  ||)  in  (17) 
is®  : 

l„(«.|sia||)=lu(f)-f:5^  (.8) 

n=l 

Thus  we  can  “subtract  out”  the  singularity  in  G{x)  by  forming 


and  then  writing  the  Green’s  function  as 

G{x)  =  G(x)  -  ^  111  jsin  ^j)  Jo  (2Ka  sin  |j)  (19') 


It  is  easy  to  see  that  G{x)  is  finite  for  all  values  of  x.  Furthermore,  both  G{x)  and 
Jo  (2«a|sin||)  in  (19')  are  periodic  and  infinitely  differentiable.  Thus  their  Fourier 
coefficients  can  be  computed  with  spectral  accuracy  using  FFT.  The  Fourier  coefficients 
of  the  Green’s  function  G(x),  gn,  will  be  computed  according  to  (19')  where  the  term 
involving  the  logarithmic  function  is  computed  by  using  convolution  sums. 


We  now  study  the  non-smoothness  of  the  normal  derivative  of  the  Green’s  function 

dG 

-^{x).  The  asymptotic  series  of  the  Bessel  functions  of  first  order  for  small  argument 
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Then 


^(i)  = '  (2,ca  |sm  -  )  sm-| 


“  ~T  l*'“  f  I)  ^  h™  f  I)]  h'“  i 
= '4^ + ^ H“  II)  •'■  (^“  1““  II)  l“”  II + 

Thus  although  ^  is  a  finite  function,  due  to  the  logarithmic  function  appearing  in 
the  second  term  shown  in  (20),  it  does  not  have  a  smooth  second  derivative  at  i  =  0. 
For  this  reason,  its  Fourier  approximation  will  converge  only  at  the  rate  of  \/N^. 

The  Fourier  coefficients  of  however,  can  be  found  easily  using  the  relation  to 
gn  given  in  Appendix  A.  In  particular,  we  have 

f -^(<7n+l -yn-l)  n/0,-f,f-l 


hn  = 


-^{92  -  90)  -  -51 
4  a 


4n  ^-T+i 


4n  T 


n-2 


n  =  0 


- 

~  2 


n  =  f -1 


Thus  it  is  only  necessary  to  compute  5n,  the  Fourier  coefficients  of  G(i). 


C.  Fast  Fourier  Transforms 

The  numerical  implementation  of  computing  5n  by  (19')  is  given  in  this  subsection. 
Let  us  introduce  Fourier  collocation  points 


For  convenience  of  discussion,  denote  the  following  Fourier  series  approximations 

N/2-1 

5(x)=  ^  (22a) 

n=-N/2 
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7V/2-1 

Jo  (2/ca  |sin  II)  =  Pnc""** 

n=-^/2 


(226) 


The  coefficients  of  these  expansions  are  computed  by  FFT  (backward  in  the  usual 
sense)  as  follows: 


7=0 

p„  =  ^  *^0  |sin  y  I) 

7=0 


e»nx7 


(22a') 

(226') 


in  which  G{xj)  is  computed  by  (19).  For  the  value  of  G(x)  at  i  =  0,  the  following 
limit,  obtained  from  (17),  can  be  used  : 


G(0)  =  -^  +  l 


In  addition,  we  denote  (18)  as 

OO 

In  jsm  ^  One”* 


—in* 


(22c) 


where  ag  =  In  and  a„  =  —  ^  ^  ,  for  n  ^  0. 

V  2  /  2  |n| 

Then,  by  (19'),  the  Fourier  coefficients  of  G(x)  is  computed  as 


9n  —  9n  ~ 


(23) 


where  ««  is  the  convolution  sum  : 


Ar/2-i 

^  ]  PmO’n—m 


(24) 


m——N/2 

We  note  that  the  convolution  sums  in  (24)  require  N  multiplications  for  each  ««. 
Thus  the  total  operations  for  the  convolution  sums  are  of  order  0{N^).  This  cost  can 
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be  reduced  considerably  to  0{N  log2  N)  by  the  use  of  a  pseudospectral  transformation 
method  with  de-aliasing  techniques^’*.  For  completeness,  evaluation  of  (24)  with  a 
“padding”  de-aliasing  technique  is  given  in  Appendix  B. 


III.  SPECTRAL  METHOD  FOR  NEUMANN  BOUNDARY  CONDITIONS 

We  now  discuss  the  Fourier  spectral  method  for  the  Boundary  Integral  Equation 

(8b)  of  the  Neumann  boundary  conditions.  Upon  substituting  the  truncated  Fourier 

series  of  the  layer  distribution  function  f{9)  into  (8b),  we  get 
.  iV/2-l  N/2-1  ,  2t  / 

f  E  E  /./  e-a..  =  E 

n=-Nl2  n=-Nl2  \  /  n=-Ar/2 


where  bn  are  the  Fourier  coefficients  of  the  specified  Neumann  boundary  condition. 

Again  the  integral  appearing  in  equation  (25)  is  directly  related  to  the  Fourier 
coefficients  of  and  It  is  easy  to  find  that  the  Fourier  coefficients  of  ^  are 
the  same  as  those  of  already  given  in  the  previous  section  as  hn-  The  apparent 
difficulty  here  is  with  the  second  normal  derivative  of  the  Green’s  function  It  can 

be  shown  that  this  function  is  strongly  singular  at  x  =  0  and,  indeed,  is  not  integrable 
in  the  ordinary  sense.  Fortunately  it  can  also  be  shown  that  the  integral  with  the 
second  normal  derivative  can  be  transformed  into  one  involving  tangential  derivatives 
with  reduced  singularity.  In  particular,  we  have** 


1  dn'dn 


1  2/  ^  ins] 


\  d  \  d 

where  and  --^77  represent  tangential  derivatives  on  the  boundary. 
aau  a  Oa' 

The  right  hand  side  of  (26)  is  now  integrable  in  the  sense  of  Cauchy  Principal  Value. 
To  show  this,  we  only  need  to  note  that  by  the  expression  of  the  Green’s  function  given 
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in  (6)  we  get 

=  “T  h‘“  f  I)  =  (^“  I™  II)  ^  l“”  I 

**£0  .  ®l\  sinx 

=  -^H\  ^  {2Ka  sin—  )  t— : — zr  (27) 

8  ^  V  21/  I  sin  f  I  ^  ^ 

Recalling  (20),  the  asymptotic  expression  of  for  small  x  is  found  as 

dG  sinx  ,  /  I  .  .  x|\  sinx 

w  “  S  ('“■  P'“ 2 1) V’“‘  ““ 2I) likfi + 

where  0(x)  denotes  smooth  terms  of  order  x  and  higher. 

The  singular  first  term  shown  above  is  integrable  in  the  sense  of  the  Cauchy  Prin¬ 
cipal  Value.  In  fact,  we  have 


Jo  sin^  2isign(n 


when  n  =  0 
)  when  n  ^  0 


Upon  substituting  x  =  6  —  6'  and  equating  the  coefficients  of  e*”®  ,  equation  (25) 
is  reduced  to 

y/n  +  /« 

in  which  we  have  used  the  fact  that,  for  a  circular  cylinder, 

n'  •  n  =  cos(0  -  o') 

The  integral  in  (30)  will  now  be  evaluated  through  the  Fourier  coefficients  of  each  term. 
For  the  first  term,  the  Fourier  coefficients  of  ^  are  obtained  from  the  relation 

dG  dG  ,  . 

^  =  E  “«»'■  (31) 

n=-Nl2 

where  gn  are  the  Fourier  coefficients  of  G(x)  by  (13). 
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The  Fourier  series  approximation  of  the  second  term  in  the  integral  of  (30)  can  also 
be  found  using  since  we  have 

N/2-1  N/2-1 

cos(i)G(x)  =  cai(i)  Y,  T.  9"'"'“ 

«=-JV/2  n=-N/2 

where 

129-1^+1  ”  =  -^/2 

K^n-i  +  ^n+i)  -N/2  +  1  <  n  <  N/2  -  2  (33) 

\9n_2  n  =  Nf2-l 

Hence  equation  (30)  is  reduced  to  the  following  algebraic  equations 

y/n  +  2jra/„  (^-^9n  +  =  K  (34) 

for  -N/2  <  n  <  7V/2  -  1. 

Therefore,  the  Fourier  coefficients  of  the  layer  distribution  function  f{6)  for  the 
Neumann  boundary  conditions  are  obtained  explicitly  as 

fn  =  ^  j  Y 

y  +  2X0  l—^9n  +  -  t»/^n  j 

where  gn  and  hn  are  computed  by  (23),  (33),  and  (21),  respectively. 

We  point  out,  however,  that  gn  as  given  by  (33)  and,  indeed,  hn  of  (21),  are  not 
exact  for  n  =  —N/2  and  N/2  —  1,  owing  to  a  truncated  series  of  G{x)  in  the  compu¬ 
tation.  Whereas  it  is  possible  to  compute  these  two  coefficients  exactly,  the  resulting 
error  in  the  last  two  coefficients  of  /«  is  negligible  because  bn,  in  the  numerator,  de¬ 
cays  exponentially  as  for  smooth  boundary  conditions.  That  is,  /«  for  n  =  —N/2  and 
N/2  —  1  are  necessarily  negligibly  small  if  N  is  sufficiently  large.  For  simplicity  and 
practicality,  (21)  and  (33)  are  retained  in  the  numerical  calculations. 
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IV.  NUMERICAL  EXAMPLES 

In  this  section,  numerical  results  of  a  plane  wave  scattering  by  a  circular  cylinder 
are  presented.  The  incident  wave  is  assumed  to  be 

A.  _  AkX 

<pi  =  c 

The  scattered  wave,  <f>,  satisfies  the  Helmholtz  equation  (1).  The  boundary  conditions 
considered  here  are  the  Dirichlet  type  <f>  —  — and  the  Neumann  type  ^ 

The  solutions  for  the  scattered  field  are  obtained  by  the  layer  representation  (2)  as 

^(r')  =  -  i<,a)  maM 


The  above  integral  can  be  easily  evaluated  directly  using  FFT,  since  the  Green’s  func¬ 
tion  has  no  singularity  for  points  lying  outside  of  the  boundary.  The  details  are  omitted 
here. 

For  plane  incident  waves,  exact  solution  is  given  by  infinite  series  of  the  Bessel  and 
Hankel  functions^.  Ovir  purpose  here  is  to  demonstrate  the  exponential  rate  of  conver¬ 
gence  of  the  numerical  solutions.  We  emphasize  again  that  the  numerical  formulation 
applies  to  any  form  of  the  incident  waves.  Due  to  its  simplicity,  a  sample  FORTRAN 
program  is  listed  in  Appendix  C. 

In  numerical  calculations,  the  radius  of  the  cylinder,  a,  is  taken  to  be  1  and  also 
1/  =  1.  Computations  for  /ca  =  1,  10  and  100  have  been  carried  out.  In  Tables  I  to 
IV,  numerical  values  of  the  layer  distribution  function  f{0)  and  the  scattered  function 
<i>  at  far  field  are  given  for  selected  points  in  space.  Exact  values  at  far  field  are  also 
shown  in  the  tables.  Clearly  as  the  number  of  Fourier  collocation  points  increases,  the 
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numerical  solution  converges  exponentially  faist.  Significant  improvements  in  accuracy 
are  observed  with  relatively  small  increase  of  the  number  of  datum  points.  This  is 
often  true  for  spectral  methods  in  general.  The  error  decreases  dramatically  when  the 
number  of  points  is  large  enough  to  resolve  the  basic  features  of  the  solution. 

The  corresponding  layer  distribution  function  f{6)  is  plotted  in  Figures  1  to  6 
for  the  Dirichlet  and  Neumann  boundary  conditions  for  /ca  =  1,  10  and  100.  These 
graphs  demonstrate  again  the  remarkable  accuracy  of  the  Fourier  spectral  methods 
with  relatively  small  number  of  datum  points. 

Fcir  field  scattered  intensities,  computed  as  lr|^^,  are  plotted  in  Figure  8  and  9  for 
the  Dirichlet  and  Neumann  boundary  conditions,  respectively. 


V.  CONCLUSIONS 

A  fast  numerical  solution  of  wave  scattering  by  a  circular  cylinder  has  been  pre¬ 
sented.  It  is  shown  that  by  properly  removing  the  non-smoothness  of  the  integral  ker¬ 
nels  of  the  Boundary  Integral  Equations,  Spectrally  accurate  numerical  solutions  are 
obtained.  The  numerical  error  decays  exponenti2dly  as  the  number  of  datum  points 
increase.  This  implies  that  the  present  method  requires  significantly  fewer  points  for 
2u:hieving  a  given  accur2M:y  in  comparison  with  previous  numerical  approaches.  The 
present  method  is  also  easy  to  implement. 

Moreover,  the  combined  single  and  double-layer  formulation  of  the  Bound2ny  Inte¬ 
gral  Equations  ensures  the  uniqueness  of  the  numerical  solution  for  all  frequencies.  It 
is  shown  that  the  hypersingularity  of  the  Boundary  Integral  Equations  can  be  handled 
easily  in  the  spectral  method. 
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APPENDIX 

A.  Exact  expressions  of  gn  and  /i„ 

In  this  appendix  we  derive  the  exact  an2dytic  expressions  for  the  Fourier  coefficients 
of  c;(x)  «.d  ^(x). 

It  can  be  shown  that,  e.g.  by  (7.2.51)  of  ref.  6, 

//J')(2Ka|sin||)=  (Al) 

m=— oo 


_ 

27r 


^0 

^"i/fW(2.o|sio||)e-x 

Moreover,  for  n  ^  0,  using  integration  by  part  and  (Al), 
—  I  -5-(x)c‘’**dx 


2ir  Jq  dn ' 


4xn 


r 


^2x0 

.  z  > 
^^°2> 

)  sm- 

^2Ka  1 

"“IDhIl 

e’"*dx 


2 


sin  ze****  dx 


=  ^/  [  £  H^m\Ka.)Jm{Ka,y 

Lm=— OO  j 

K^a 

-  ”‘^«n+l  ~9n-l) 

where  use  has  been  made  of  the  formula^ 

^  [xAf<"(x)]  =  xff<'>(x) 
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For  n  =  0,  further  calculations  show  that 


2  ’ 

ho  =  [//5*^(ica)J2(ica)  -  //J‘\Ka)Jo('ta)]  - 


4a 


K^a  1 

=  -7“ W  “  SO) - 91 

4  a 


B.  Evaluation  of  convolution  sums 


An  algorithm  of  computing  convolution  sums  Un  with  0{N  log2  N)  operations  is 
shown  below®. 


Let  M  >  3N  and 


(j  =  2irj/M  ,  j  =  0, 1, 2,  ...Af  -  1 


Compute  the  following  using  FFT  for  j  =  0, 1,2,  ...Af  —  1  : 

lli/2-l 

m=-Mf2 

Uf2-l 


where 


ms:— hi  12 


{Om  —N  <  m  <  N  —  \ 
0  other 


and  form  the  product 


Pm  -N/2  <m<Nl2-\ 


other 


Uj  =  AjPj 


Then  the  convolution  sum  Un  is  the  (backward)  FFT  of  Uj  as  follows 


,  M-i 
)=0 
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for  -JV/2  <  n  <  iV/2  -  1. 


C.  FORTRAN  program 

A  FORTRAN  program  of  implementing  the  Fourier  spectral  method  is  listed  below. 
(The  routines  cftti,  cfttf  and  cfttb  denote  initialization,  forward  and  backward  FFT 
transforms  respectively.) 

program  circla 


c  n  :  nnabar  of  points;  isoftsl  :  Dirichlat  B.C.;  isoltsO  :  lanmann  B.C. 

c************************************************************************* 
paraaatar (n>32 , ak»10 . 0 , isof t*! ,ota«l . O.nlen-l .nlialf sB/2,m«3*B, 

>  maf lost (n) .pis3 . 141692«S3S8979324,aalar«0 . 57721666490153286) 
coaplaz  b(0:al) ,fn(0:al) ,gbar(0:nl) ,gn(0:Bl) ,bn(0:Bl) ,p(0:Bl) , 

>  gtildo(0:nl),aB(0:B-l),pa(0:n-l),*saTo(2000),ssaTa2(2000).ai.phi 
oi3(0, 0.1.0) 

call  clftKn.vsava) 
call  gatbc(n,ak,b,ai,pi) 
call  cfftf (n.b.asava) 
do  10  j*0,n-l 

tBp«2 . 0«ak*abs(sln(pi*tloat  ( j  )/m)  ) 

If(j.aq.O)  than 
gbar (0 ) a-anlar/2 . O/pi+ai/4 . 0 
p(0)ol.0 

also 

gbar(  j  )3ai/4. 0«(bas  j0(tBp)'»’ai*basp0(ti^)  ) 

>  -»'0.5*alog(tap/2.0)abasj0(ti9)/pi 
p(j)=basj0(tnp) 

andif 

10  continna 

call  cfftb(n,gbar,«sava) 
call  cfftb(n,p,«sava) 
aB(0)salog(ak/2 . 0) 
aB(2*n)  =-l .  0/2 . 0/m 
do  21  isl,n-l 

amC i ) s- 1 . 0/2 . 0/float ( i ) 

21  aB(2*n+i)=1.0/2.0/float(i-n) 
do  22  isO.nhalf-l 

Fai(i)»pU) 

22  pm(6*nbalf+i)=p(nhalf-*-i) 
call  cffti(B,asaTa2) 

call  cfftf (B,aB,«saTo2) 
call  cfftf (B,pm,vsaTo2) 
do  23  j=0,B-l 

23  pai(j)aaB(j)*pB(j) 
call  cfftb(B.pa,«sava2) 
do  31  isO.nhalf-l 
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gB(i)«gbar(i)-0 . 6*|a(i)/f lo&t(B)/pi 

31  gB(ali»lf'«'i)agb«r(iihalf'<-l)-0.6*pa(5*iilMlf*i)/float(B)/pl 
ha(0)-ak**2/4. 0*(gB(2)-ga(0> )-gR(l) 

hn(iiK«lf-l )  ■•k**2/4 . 0/float  (ahalf -1  )*ga(alialf-2) 
ha(Bhalf  )  «a]co*2/4 . 0/float  (ahalf  )  oga  (ahalf  4^  1 ) 
ha(a-l  )*>ak*o2/4 . 0*(ga(0)-ga(a-2)  ) 
gtildo(0)«0 . 5*(ga( 1 )4ga(a-l ) ) 
gtildo(ahalf -1 ) *0 . S*ga(ahalf-2) 
gtilda(ahalf )«0.S*ga(ahalf4l) 
gtilda(a-l)>0.S*(ga(0)4ga(a-2)) 
do  32  l>l,a-2 
itno«i 

if (i. go. ahalf )  ltno«i-a 

if (i.oq.ahalf-l.or.i.oq.ahalf)  go  to  32 

ha( i ) •-ak*02/4 . 0/float ( itrno )• (ga( i4l ) -ga ( i- 1 ) ) 

gtildo(i)«0.6*(ga(i  ■l)4ga(i4l)) 

32  coatiaua 

do  40  i30,a-l 

if (isoft.oq. 1)  thoa 

f a ( i ) >b( i ) / ( 0 . 5*xn42 . 0*pi* (ha( i ) -oi*ota*ga ( i ) ) ) 

olao 

itrao^i 

if (i. go. ahalf )  itrao«i-a 

fa(i)«b(i)/(0. 6*oi*ata*ra+2 . 0*pi*(-f loat (itrao) **2*ga( i ) 
>  4ak**2*gtildo(i)-oi*ota*ha(i))) 

oadif 

40  coatiaua 

c*********************************************************** 
c  Tho  follouiag  is  to  fiad  phi  at  far  fiold  r»rO 
c*********************************************************** 

rOslO.O 

apoiata4 

do  70  ii»l,Bpoiat 

0 j-2 . 0*pi*f loat ( ii-1 ) /float (apoiat) 
do  71  j»0,a“l 

thotas2 . 0*pi*f loat ( j ) /ra 
r j  ssqrt ( 1 . 0+r0*r0-2 . 0*r0*co« (thata-s j ) ) 
djsl .0-r0*cos(thota-sj) 
tap=ak*rj 

gB( j ) «oi/4 . 0* (boa j  0 (tiq>)4ei*boay0 (tap) ) 

71  hB( j ) s-oi*ak/4 . 0*(boa j 1 (tap)4oi*beayl (tap) ) *dj/r j 
call  cfftb(a,ga,waaTo) 

call  cfftb(a,ha,«aaTo) 
phi=0 . 0 
do  72  i«0,a-l 

72  phisphi42.0*pi*fa(i)*(hn(i)-oi*ota*go(i))/ra 

70  «rito(3,100)  rO,aj ,phi,caba(phi) 

100  foraat('  r0>'.ol5.6.'  thota»' ,ol5.6/'  pbi=’ .3ol7. 10) 

999  atop 
oad 
c 

aubrontiao  gotbc(a,ak,b,oi,pi) 
coaploz  oi,b(0:  a-l),tap 
do  10  j=0,a-l 
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tap3ai*ak*cos (2 . 0*pi*l loat ( j )/f loat (n) ) 
10  b(j)=-c«xp(tmp) 
return 
end 
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TABLE  I 

Values  of  the  layer  distribution  function  f{9)  at  selected  points 
on  the  boundary.  Dirichlet  boundary  condition. 


N 

^  =  0" 

5  =  90* 

e 

O 

00 

II 

Error  | 

1 _ ''iri _  1 

4 

1.101447573 

1.102982967 

1.124378820 

10-2 

8 

1.113205176 

1.095419894 

1.146430615 

10-2 

16 

1.112753432 

1.094877536 

1.145739275 

io-« 

24 

1.112753420 

1.094877525 

1.145739263 

10->2 

Ka  =  10 

24 

4.590213453 

6.904710445 

5.180354736 

10-2 

32 

4.546357630 

6.901732036 

5.132718905 

10-2 

48 

4.545461066 

6.901500667 

5.132515158 

10-* 

56 

4.545461055 

6.901500659 

5.132515156 

10-^2 

Ka  =  100 

224 

20.64255659 

6.841653547 

18.93255934 

10-2 

256 

20.64325731 

6.842244857 

18.93221646 

10-^ 

512 

20.64325733 

6.842244863 

18.93221644 

10-»2 
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TABLE  II 

Values  of  the  scattered  function  (j>  at  selected  points 
at  far  field  r  =  10a.  Dirichlet  boundary  condition. 


N 

II 

o 

o 

II 

CO 

o 

o 

e  =  180“ 

Error 

Ka  =  1 

4 

0.4146449903 

0.2787718545 

0.1852248716 

10-2 

8 

0.4224209076 

0.2612785029 

0.2551151985 

lo-** 

16 

0.4224153154 

0.2613031445 

0.2552183381 

10-10 

Exact 

0.4224153154 

0.2613031445 

0.2552183381 

Ka  =  10 

24 

0.8255952003 

0.1969679200 

0.1864749710 

10-2 

32 

0.8285176644 

0.1953580665 

0.2300067055 

10-'* 

48 

0.8285110664 

0.1953543814 

0.2300939707 

10-10 

Exact 

0.8285110664 

0.1953543814 

0.2300939707 

Ka  =  100 

224 

0.8562228283 

0.1881301853 

0.2295232548 

10"^ 

256 

0.8562289911 

0.1881326409 

0.2294229274 

10-10 

Exact 

0.8562289911 

0.1881326409 

0.2294229274 
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TABLE  III 

Values  of  the  layer  distribution  function  f{0)  at  selected  points 
on  the  boundary.  Neumann  boundary  condition. 


N 

II 

o 

o 

II 

CO 

o 

o 

II 

00 

o 

o 

Error 

/ca  =  1 

4 

1.035182633 

0.3028073027 

0.8616587030 

10-^ 

8 

1.200134116 

0.3972281648 

0.8518411247 

10-2 

16 

1.199187560 

0.3963806796 

0.8495643896 

lo-"^ 

24 

1.199187560 

0.3963806589 

0.8495643587 

10-12 

Ka  =  10 

24 

0.6004486353 

0.4814454225 

1.362228889 

10-1 

32 

0.6274625969 

0.6575899642 

1.577833267 

10-2 

48 

0.6302381163 

0.6567081358 

1.460119442 

10-^ 

56 

0.6302381517 

0.6567081198 

1.460119455 

10-12 

Ka  =  100 

224 

0.2185547081 

1.282490390 

2.054272775 

10-2 

256 

0.2157948725 

1.283008634 

10-^ 

512 

0.2157947803 

1.283008643 

2.057913072 

10-12 

24 


TABLE  IV 

Values  of  the  scattered  function  <i>  at  selected  points 
at  far  field  r  =  10a.  Neumann  boundary  condition. 


N 

0  =  0" 

0  =  90" 

O 

O 

00 

II 

Error 

'  /ca  =  1 

4 

0.1583300606 

0.1690204144 

0.1619964200 

10"^ 

8 

0.1732916160 

0.1563414831 

0.2312523394 

lO"'* 

16 

0.1733358919 

0.1563260243 

0.2313583724 

10-10 

Exact 

0.1733358919 

0.1563260243 

0.2313583724 

/ca  =  10 

24 

0.7679584467 

0.2167643382 

0.1574136977 

10"^ 

32 

0.7740714632 

0.1956069424 

0.2282238894 

10"^ 

48 

0.7740874173 

0.1955960691 

0.2283394143 

10-10 

Exact 

0.7740874173 

0.1955960691 

0.2283394143 

o 

o 

II 

<3 

224 

0.7688015277 

0.1871656432 

0.2295250315 

10"® 

256 

0.7688018590 

0.1871717295 

0.2293995512 

10-10 

Exact 

0.7688018590 

0.1871717295 

0.2293995512 
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Figure  1.  Layer  distribution  function  f{B)  for  na  =  \.  Dirichlet  boundary  condition. 
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Figure  2.  Layer  distribution  function  f{6)  for  /co  =  10.  Dirichlet  boundary  condition. 
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Figure  3.  Layer  distribution  function  f{d)  for  «a  =  100.  Dirichlet  boundary  condition. 
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Figure  5.  Layer  distribution  function  f{9)  for  /to  =  10.  Neumann  boundary  condition. 
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Figure  6.  Layer  distribution  function  f{0)  for  kg  =  100.  Neumann  boundary  condition. 
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Figure  8.  Directivities  of  the  far  filed  scattered  function,  Neumann  boundary  condition. 
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